R Vignette
Focusing the GWAS Lens on days to flower using latent variable phenotypes derived from global multi-environment trials
Introduction
This vignette contains the R code and analysis done for the paper:
Derek Wright, Sandesh Neupane, Jakob Butler, Raul Martinez, Jim Weller, Kirstin E Bett (2020) Focusing the GWAS Lens on days to flower using latent variable phenotypes derived from global multi-environment trials .
which is follow-up to:
Derek Wright, Sandesh Neupane, Taryn Heidecker, Teketel Haile, Clarice Coyne, Sripada Udupa, Eleonora Barilli, Diego Rubiales, Tania Gioia, Reena Mehra, Ashutosh Sarker, Rajeev Dhakal, Babul Anwar, Debashish Sarker, Albert Vandenberg, and Kirstin E. Bett. (2020) Understanding photothermal interactions can help expand production range and increase genetic diversity of lentil (Lens culinaris Medik.). Plants, People, Planet. 00:1-11.
This work done as part of the AGILE project at the University of Saskatchewan.
GWAS
# Genotypes
myG <- read.csv("myG_LDP.csv", header = F)
# Phenotypes
myY <- read.csv("myY.csv")
# CoVariates
myCV <- myY[,c("Name","b","c")]# Load library
library(GAPIT3) # devtools::install_github("jiabowang/GAPIT3",force=TRUE)
#
setwd("Results")
myGAPIT <- GAPIT(
Y = myY,
G = myG,
model = c("MLM","MLMM","FarmCPU","Blink"),
PCA.total = 4
)
# GWAS with b covariate (Results in `Results_b/`)
setwd("Results_b")
myGAPIT <- GAPIT(
Y = myY%>%select(-PTModel_b),
G = myG,
CV = myCV[,c("Name","PTModel_b")]
model = c("MLM","MLMM","FarmCPU","Blink"),
PCA.total = 0
)
# GWAS with c covariate (Results in `Results_c/`)
setwd("Results_c")
myGAPIT <- GAPIT(
Y = myY%>%select(-PTModel_c),
G = myG,
CV = myCV[,c("Name","PTModel_c")]
model = c("MLM","MLMM","FarmCPU","Blink"),
PCA.total = 0
)Prepare Post GWAS Data
# Load libraries
library(tidyverse)
library(ggpubr)
library(ggbeeswarm)
library(ggtext)
# Genotype and Phenotype
myG <- read.csv("myG_LDP.csv", header = T)
myY <- read.csv("myY.csv")
# Genotype metadata - Cluster groups for the LDP (Wright et al., 2020)
myLDP <- read.csv("myLDP.csv") %>%
mutate(Cluster = factor(Cluster))
# List of flowering time genes
myFTGenes <- read.csv("Lentil_FT_Genes.csv") %>%
rename(Chromosome=Chr) %>%
mutate(Expt = "FT genes", Model = "MLM", Position = as.numeric(Start))
# Nucleotide symbols
#myN <- read.csv("IUPAC_Nucleotide_Code.csv")
# ggplot theme
theme_AGL <- theme_bw() +
theme(strip.background = element_rect(colour = "black", fill = NA, size = 0.5),
panel.background = element_rect(colour = "black", fill = NA, size = 0.5),
panel.border = element_rect(colour = "black", size = 0.5),
panel.grid = element_line(color = alpha("black", 0.1), size = 0.5),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank())
#
threshold <- -log10(0.05/nrow(myG))
threshold2 <- -log10(0.000005)
myModels <- c("MLM", "MLMM", "FarmCPU", "Blink")
myColors_Model <- c("darkgreen","darkred","darkblue","darkslategray4")
myColors_Macro <- c("azure4","darkgreen","darkorange3","darkblue")
myColors_Cluster <- c("darkred", "darkorange3", "darkgoldenrod2", "deeppink3",
"steelblue", "darkorchid4", "cornsilk4", "darkgreen")
myMs <- c(
"Lcu.2RBY.Chr2p42556949", #"Lcu.2RBY.Chr2p42543877",
"Lcu.2RBY.Chr5p1069654", #"Lcu.2RBY.Chr5p1063138",
"Lcu.2RBY.Chr6p2528817",
"Lcu.2RBY.Chr6p307256203")#"Lcu.2RBY.Chr6p306914970") #
myGM <- read.csv("Results/GAPIT.MLM.Ro17_DTF.GWAS.Results.csv") %>%
filter(SNP %in% myMs)
#
myLabels <- c(
"PC1", "PC2", "PC3",
"*a*", "*b*", "*c*",
"Ro16_DTF", "Ro17_DTF", "Su16_DTF", "Su17_DTF", "Su18_DTF", "Us18_DTF",
"In16_DTF", "In17_DTF", "Ba16_DTF", "Ba17_DTF", "Ne16_DTF", "Ne17_DTF",
"Mo16_DTF", "Mo17_DTF", "Sp16_DTF", "Sp17_DTF", "It16_DTF", "It17_DTF",
"Su17_*Tf*", "Ba17_*Tf*", "It17_*Tf*",
"Su17_*Tb*", "Ba17_*Tb*", "It17_*Tb*",
"Su17_*Pf*", "Ba17_*Pf*", "It17_*Pf*",
"Su17_*Pc*", "Ba17_*Pc*", "It17_*Pc*")
# Figure 01
myTraits <- c(
"PCA_PC1", "PCA_PC2", "PCA_PC3",
"PTModel_a", "PTModel_b", "PTModel_c",
"Ro16_DTF", "Ro17_DTF", "Su16_DTF", "Su17_DTF", "Su18_DTF", "Us18_DTF",
"In16_DTF", "In17_DTF", "Ba16_DTF", "Ba17_DTF", "Ne16_DTF", "Ne17_DTF",
"Mo16_DTF", "Mo17_DTF", "Sp16_DTF", "Sp17_DTF", "It16_DTF", "It17_DTF",
"Su17_Tf", "Ba17_Tf", "It17_Tf",
"Su17_Tb", "Ba17_Tb", "It17_Tb",
"Su17_Pf", "Ba17_Pf", "It17_Pf",
"Su17_Pc", "Ba17_Pc", "It17_Pc")Supplemental Table 1
dropNAcol <- function (x) { x[, colSums(is.na(x)) < nrow(x)] }
GWAS_PeakTable <- function(folder = NULL, file = NULL,
g.range = 3000000, rowread = 2000) {
#
trait <- substr(file, gregexpr("GAPIT.", file)[[1]][1]+6,
gregexpr(".GWAS.Results.csv", file)[[1]][1]-1 )
model <- substr(trait, 1, gregexpr("\\.", trait)[[1]][1]-1 )
trait <- substr(trait, gregexpr("\\.", trait)[[1]][1]+1, nchar(trait) )
expt <- substr(trait, 1, gregexpr("_", trait)[[1]][1]-1)
trait <- substr(trait, gregexpr("_", trait)[[1]][1]+1, nchar(trait) )
#
output <- NULL
#
rr <- read.csv(paste0(folder, file), nrows = rowread) %>%
mutate(`-log10(p)` = -log10(P.value))
#
if(sum(colnames(rr)=="nobs")>0) { rr <- select(rr, -nobs) }
#
if(!is.null(threshold2)) {
rx <- rr %>% filter(-log10(P.value) >= threshold2)
} else{ rx <- rr %>% filter(-log10(P.value) > threshold) }
#
if(nrow(rx) == 0) {
rx[1,] <- NA
rx[1,"Chromosome"] <- 1
output <- bind_rows(output, rx)
} else {
while(nrow(rx) > 0) {
rp <- rx %>% group_by(Chromosome) %>%
top_n(., n = 1, `-log10(p)`) %>% ungroup()
output <- bind_rows(output, rp)
# i<-2
for(i in 1:nrow(rp)) {
g.range1 <- rp$Position[i] - g.range
g.range2 <- rp$Position[i] + g.range
if(g.range1 < 0) { g.range1 <- 0 }
rx <- rx[rx$Chromosome != rp$Chromosome[i] |
rx$Position < g.range1 |
rx$Position > g.range2,]
}
if(!is.null(threshold2)) {
rx <- rx %>% filter(-log10(P.value) > threshold2)
} else{ rx <- rx %>% filter(-log10(P.value) > threshold) }
}
}
output <- output %>%
mutate(Model = model, Trait = trait, Expt = expt,
FT_Genes = NA, FT_Pos = NA, FT_End = NA,
FT_Distance = NA, FT_Distances = NA)
#
if(sum(!is.na(output$Position)) > 0) {
for(i in 1:nrow(output)) {
g.range1 <- output$Position[i] - g.range
g.range2 <- output$Position[i] + g.range
ftg <- myFTGenes[myFTGenes$Chromosome == output$Chromosome[i] &
myFTGenes$Position > g.range1 &
myFTGenes$Position < g.range2,]
ftg <- ftg %>%
mutate(Distance = abs(Position - output$Position[i])) %>%
arrange(Distance)
output$FT_Genes[i] <- paste(ftg$Name, collapse = " ; ")
output$FT_Pos[i] <- paste(ftg$Position, collapse = " ; ")
output$FT_End[i] <- paste(ftg$End, collapse = " ; ")
output$FT_Distances[i] <- paste(ftg$Distance, collapse = " ; ")
output$FT_Distance[i] <- ftg$Distance[1]
}
}
output %>%
select(SNP, Chromosome, Position, Model, Expt, Trait,
P.value, `-log10(p)`, effect, MAF=maf,
Rsquare.of.Model.without.SNP, Rsquare.of.Model.with.SNP,
FDR_Adjusted_P.values,
FT_Genes, FT_Pos, FT_End, FT_Distance, FT_Distances) %>%
filter(!(`-log10(p)` < threshold & Model != "MLM"))
}
#
myP <- NULL
#
fnames <- grep(".GWAS.Results", list.files("Results"))
fnames <- list.files("Results")[fnames]
for(i in fnames) {
myP <- bind_rows(myP, GWAS_PeakTable(folder = "Results/", file = i) %>% dropNAcol())
}
#
fnames <- grep(".GWAS.Results", list.files("Results_b"))
fnames <- list.files("Results_b")[fnames]
for(i in fnames) {
myP <- bind_rows(myP, GWAS_PeakTable(folder = "Results_b/", file = i) %>%
mutate(Trait = paste0(Trait,"-b")) %>% dropNAcol())
}
#
fnames <- grep(".GWAS.Results", list.files("Results_c"))
fnames <- list.files("Results_c")[fnames]
for(i in fnames) {
myP <- bind_rows(myP, GWAS_PeakTable(folder = "Results_c/", file = i) %>%
mutate(Trait = paste0(Trait,"-c")) %>% dropNAcol())
}
#
myP <- myP %>% filter(!is.na(SNP)) %>%
arrange(Chromosome, Position, P.value, Trait) %>%
mutate(Chromosome = factor(Chromosome, levels = 1:7),
Model = factor(Model, levels = myModels))
#
write.csv(myP, "Supplemental_Table_01.csv", row.names = F)Figure 1
gg_gwas_summary <- function(traits, labels = traits,
hlines, width = 12, height = 8,
title = NULL, filename) {
#
me <- c("<b style='color:black'>{ExptTrait}</b>",
"<b style='color:darkgreen'>{ExptTrait}</b>",
"<b style='color:darkorange3'>{ExptTrait}</b>",
"<b style='color:darkblue'>{ExptTrait}</b>")
#
e1 <- c("Ro16","Ro17","Su16","Su17","Su18","Us18")
e2 <- c("Mo16","Mo17","Sp16","Sp17","It16","It17")
e3 <- c("In16","In17","Ba16","Ba17","Ne16","Ne17")
myMacroEnvs <- c("Multi Environment","Temperate", "South Asia", "Mediterranean")
#
xx <- myP %>%
filter(paste(Expt, Trait, sep = "_") %in% traits) %>%
mutate(ExptTrait = paste(Expt, Trait, sep = "_"),
ExptTrait = plyr::mapvalues(ExptTrait, traits, labels),
Model = factor(Model, levels = myModels),
MacroEnv = NA,
MacroEnv = ifelse(Expt %in% e1, "Temperate", MacroEnv),
MacroEnv = ifelse(Expt %in% e2, "Mediterranean", MacroEnv),
MacroEnv = ifelse(Expt %in% e3, "South Asia", MacroEnv),
MacroEnv = ifelse(is.na(MacroEnv), "Multi Environment", MacroEnv),
cat_col = NA,
cat_col = ifelse(MacroEnv == "Multi Environment", glue::glue(me[1]), cat_col),
cat_col = ifelse(MacroEnv == "Temperate", glue::glue(me[2]), cat_col),
cat_col = ifelse(MacroEnv == "South Asia", glue::glue(me[3]), cat_col),
cat_col = ifelse(MacroEnv == "Mediterranean", glue::glue(me[4]), cat_col),
cat_col = factor(cat_col, levels = rev(cat_col[match(labels, ExptTrait)])),
ExptTrait = factor(ExptTrait, levels = rev(labels)),
MacroEnv = factor(MacroEnv, levels = myMacroEnvs)) %>%
arrange(ExptTrait)
x1 <- xx %>% filter(`-log10(p)` > threshold)
x2 <- xx %>% filter(`-log10(p)` < threshold)
#
mp <- ggplot(x1, aes(x = Position / 100000000, y = cat_col,
shape = Model, fill = MacroEnv)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_point(data = x2, size = 1, color = "black", alpha = 0.5,
aes(key1 = SNP, key2 = FT_Genes, key3 = FT_Distances)) +
geom_point(size = 2, color = "black", alpha = 0.5,
aes(key1 = SNP, key2 = FT_Genes, key3 = FT_Distances)) +
geom_hline(yintercept = hlines, alpha = 0.7) +
facet_grid(. ~ Chromosome, drop = F, scales = "free_x", space = "free_x") +
scale_fill_manual(values = myColors_Macro, guide = "none") +
scale_shape_manual(values = c(22:25), breaks = myModels) +
scale_y_discrete(drop = F) +
scale_x_continuous(breaks = 0:6, minor_breaks = 0:6) +
theme_AGL +
theme(legend.position = "bottom",
axis.text.y = element_markdown(),
strip.text.y = element_markdown(),
plot.title = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 4, fill = "black"))) +
labs(title = title, y = NULL, x = "100 Mbp")
#
ggsave(paste0(filename, ".png"), mp, width = width, height = height)
#
mp <- plotly::ggplotly(mp)
htmlwidgets::saveWidget(plotly::as_widget(mp),
paste0(filename,"_plotly.html"),
knitrOptions = list(fig.width = width, fig.height = height),
selfcontained = T)
}
#
gg_gwas_summary(traits = myTraits, labels = myLabels,
hlines = c(3.5,6.5,9.5,12.5,18.5,24.5,30.5,33.5),
width = 12, height = 9, filename = "Figure_01")
#
gg_gwas_summary(title = "Covariate = *b*",
traits = paste0(myTraits, "-b"),
labels = myLabels,
hlines = c(3.5,6.5,9.5,12.5,18.5,24.5,30.5,33.5),
width = 12, height = 9, filename = "Additional/Additional_Figure_01")
#
gg_gwas_summary(title = "Covariatie = *c*",
traits = paste0(myTraits, "-c"),
labels = myLabels,
hlines = c(3.5,6.5,9.5,12.5,18.5,24.5,30.5,33.5),
width = 12, height = 9, filename = "Additional/Additional_Figure_02")Figure 2
expttraits <- c("Su17_DTF", "Ba17_DTF", "Ne17_DTF", "It17_DTF")
expts <- expttraits
for(i in 1:length(expts)) {
expts[i] <- substr(expts[i], 1, gregexpr("_", expts)[[i]][1]-1)
}
#
yy <- myY %>%
gather(ExptTrait, Value, 2:ncol(.)) %>%
filter(ExptTrait %in% expttraits) %>%
mutate(Model = factor("Days To Flower"),
ExptTrait = factor(ExptTrait, levels = expttraits),
Expt = plyr::mapvalues(ExptTrait, expttraits, expts),
Expt = factor(Expt, levels = expts))
#
xx <- NULL
for(i in expttraits) {
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results/"))
fnames <- list.files("Results/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results/", j)) %>%
mutate(Model = mod, ExptTrait = i,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xx <- bind_rows(xx, xj)
}
}
#
xx <- xx %>%
mutate(Model = factor(Model, levels = myModels),
ExptTrait = factor(ExptTrait, levels = expttraits),
Expt = plyr::mapvalues(ExptTrait, expttraits, expts),
Expt = factor(Expt, levels = expts)) %>%
arrange(desc(Model))
#
xx <- xx %>% mutate(`-log10(p)` = ifelse(`-log10(p)` > 10, 10, `-log10(p)`))
#
x1 <- xx %>% filter(`-log10(p)` < threshold2)
x2 <- xx %>% filter(`-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(`-log10(p)` > threshold)
#
mp1 <- ggplot(yy, aes(x = Value)) +
geom_histogram(fill = "darkgreen", binwidth = 2) +
facet_grid(Expt ~ "DTF") +
theme_AGL +
xlim(c(35,160)) +
labs(x = "Days after planting", y = "Count")
#
mybreaks <- c(0, 3, threshold2, threshold, 8, 10)
mylabels <- c(0, 3, round(threshold2,1), round(threshold,1), 8, ">10")
#
mp2 <- ggplot(x1, aes(x = Position / 100000000, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
facet_grid(Expt ~ Chromosome, scales = "free", space = "free_x") +
scale_x_continuous(breaks = 0:7) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(legend.position = "none",
strip.text.y = element_blank() ,
strip.background.y = element_blank(),
axis.title.y = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "100 Mbp")
#
mp3 <- ggplot(x1, aes(x = `-log10(p)_Exp`, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_point(size = 0.5, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25, color = alpha("white", 0)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_abline() +
scale_x_continuous(breaks = mybreaks, labels = mylabels) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
facet_grid(Expt ~ "QQ", scales = "free_y") +
theme_AGL +
theme(legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank()) +
labs(y = NULL, x = "Expected")
#
mpl <- get_legend(mp2, position = "bottom")
mp <- ggarrange(mp1, mp2, mp3, ncol = 3, align = "h",
widths = c(2,7, 1.5), labels = c("a)", "b)", NULL),
legend.grob = mpl, legend = "bottom", common.legend = T)
#
ggsave("Figure_02.png", mp, width = 12, height = 7)Figure 3
expttraits <- c("PCA_PC1", "PTModel_b", "PCA_PC2", "PTModel_c", "PCA_PC3", "Su17_Tb")
expts <- c("PC1", "*b*", "PC2", "*c*", "PC3", "Su17_*Tb*")
#
yy <- myY
for(i in expttraits) { yy[,i] <- scales::rescale(yy[,i], c(0, 1)) }
yy <- yy %>%
gather(ExptTrait, Value, 2:ncol(.)) %>%
filter(ExptTrait %in% expttraits) %>%
mutate(ExptTrait = plyr::mapvalues(ExptTrait, expttraits, expts),
ExptTrait = factor(ExptTrait, levels = expts))
#
xx <- NULL
for(i in expttraits) {
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results/"))
fnames <- list.files("Results/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results/", j)) %>%
mutate(Model = mod, ExptTrait = i,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xx <- bind_rows(xx, xj)
}
}
#
xx <- xx %>%
mutate(Model = factor(Model, levels = myModels),
ExptTrait = plyr::mapvalues(ExptTrait, expttraits, expts),
ExptTrait = factor(ExptTrait, levels = expts)) %>%
arrange(desc(Model))
#
xx <- xx %>% mutate(`-log10(p)` = ifelse(`-log10(p)` > 10, 10, `-log10(p)`))
#
x1 <- xx %>% filter(`-log10(p)` < threshold2)
x2 <- xx %>% filter(`-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(`-log10(p)` > threshold)
#
mp1 <- ggplot(yy, aes(x = Value)) +
geom_histogram(fill = "darkgreen") +
facet_grid(ExptTrait ~ "Phenotype", scales = "free") +
theme_AGL +
theme(strip.text.y = ggtext::element_markdown()) +
labs(x = "Scaled Value", y = "Count")
#
mybreaks <- c(0, 3, threshold2, threshold, 8, 10)
mylabels <- c(0, 3, round(threshold2,1), round(threshold,1), 8, ">10")
#
mp2 <- ggplot(x1, aes(x = Position / 100000000, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
facet_grid(ExptTrait ~ Chromosome, scales = "free", space = "free_x") +
scale_x_continuous(breaks = 0:7) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(legend.position = "none",
strip.text.y = element_blank() ,
strip.background.y = element_blank(),
axis.title.y = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "100 Mbp")
#
mp3 <- ggplot(x1, aes(x = `-log10(p)_Exp`, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_point(size = 0.5, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25, color = alpha("white", 0)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_abline() +
scale_x_continuous(breaks = mybreaks, labels = mylabels) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
facet_grid(ExptTrait ~ "QQ", scales = "free_y") +
theme_AGL +
theme(legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
strip.text.y = ggtext::element_markdown()) +
labs(y = NULL, x = "Expected")
#
mpl <- get_legend(mp2, position = "bottom")
mp <- ggarrange(mp1, mp2, mp3, ncol = 3, widths = c(2,7,1.5),
align = "h", labels = c("a)", "b)", NULL),
legend.grob = mpl, legend = "bottom", common.legend = T)
#
ggsave("Figure_03.png", mp, width = 12, height = 10)Figure 3 Subplots
fig3_Subplot <- function(expttraits = c("PC1", "*b*")) {
#
x1 <- xx %>% filter(ExptTrait %in% expttraits, `-log10(p)` < threshold2)
x2 <- xx %>% filter(ExptTrait %in% expttraits, `-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(ExptTrait %in% expttraits, `-log10(p)` > threshold)
#
mp1 <- ggplot(yy %>% filter(ExptTrait %in% expttraits), aes(x = Value)) +
geom_histogram(fill = "darkgreen") +
facet_grid(ExptTrait ~ "Phenotype", scales = "free") +
theme_AGL +
theme(strip.text.y = ggtext::element_markdown()) +
labs(x = "Scaled Value", y = "Count")
#
mybreaks <- c(0, 3, threshold2, threshold, 8, 10)
mylabels <- c(0, 3, round(threshold2,1), round(threshold,1), 8, ">10")
#
mp2 <- ggplot(x1, aes(x = Position / 100000000, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
facet_grid(ExptTrait ~ Chromosome, scales = "free", space = "free_x") +
scale_x_continuous(breaks = 1:7) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(legend.position = "none",
strip.text.y = element_blank() ,
strip.background.y = element_blank(),
axis.title.y = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "100 Mbp")
#
mp3 <- ggplot(x1, aes(x = `-log10(p)_Exp`, y = `-log10(p)`, shape = Model,
fill = Model)) +
geom_point(size = 0.5, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25, color = alpha("white", 0)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_abline() +
scale_x_continuous(breaks = mybreaks, labels = mylabels) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
facet_grid(ExptTrait ~ "QQ", scales = "free_y") +
theme_AGL +
theme(legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
strip.text.y = ggtext::element_markdown()) +
labs(y = NULL, x = "Expected")
#
mpl <- get_legend(mp2, position = "bottom")
mp <- ggarrange(mp1, mp2, mp3, ncol = 3, widths = c(2,7,1.5),
align = "h", labels = c("a)", "b)", NULL),
legend.grob = mpl, legend = "bottom", common.legend = T)
#
mp
}
ggsave("Additional/Figure_03_1.png",
fig3_Subplot(expttraits = c("PC1", "*b*")),
width = 12, height = 4)
ggsave("Additional/Figure_03_2.png",
fig3_Subplot(expttraits = c("PC2", "*c*")),
width = 12, height = 4)
ggsave("Additional/Figure_03_3.png",
fig3_Subplot(expttraits = c("PC3", "Su17_*Tb*")),
width = 12, height = 4)Figure 4
expttraits <- myTraits
expts <- expttraits
for(i in 1:length(expts)) {
expts[i] <- substr(expts[i], 1, gregexpr("_", expts)[[i]][1]-1)
}
#
# i<-"Ro17_DTF"
#expttraits <- c("Ro17_DTF","Sp17_DTF")
for(i in expttraits) {
x1 <- NULL; xb <- NULL; xc <- NULL
#
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results/"))
fnames <- list.files("Results/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results/", j)) %>%
mutate(Model = mod, Facet = i,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
x1 <- bind_rows(x1, xj)
}
#
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results_b/"))
fnames <- list.files("Results_b/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results_b/", j)) %>%
mutate(Model = mod, Facet = "CV = *b*",
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xb <- bind_rows(xb, xj)
}
#
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results_c/"))
fnames <- list.files("Results_c/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results_c/", j)) %>%
mutate(Model = mod, Facet = "CV = *c*",
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xc <- bind_rows(xc, xj)
}
#
xx <- bind_rows(x1, xb, xc) %>%
mutate(Model = factor(Model, levels = myModels),
Facet = factor(Facet, levels = c("CV = *c*", i, "CV = *b*"))) %>%
arrange(Facet, desc(Model))
#
xx <- xx %>% mutate(`-log10(p)` = ifelse(`-log10(p)` > 10, 10, `-log10(p)`))
#
x1 <- xx %>% filter(`-log10(p)` < threshold2)
x2 <- xx %>% filter(`-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(`-log10(p)` > threshold)
#
mybreaks <- c(0, 3, threshold2, threshold, 8, 10)
mylabels <- c(0, 3, round(threshold2,1), round(threshold,1), 8, ">10")
#
mp1 <- ggplot(x1, aes(x = Position / 100000000, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
facet_grid(Facet ~ Chromosome, scales = "free", space = "free_x") +
scale_x_continuous(breaks = 0:7) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(strip.text.y = element_blank() ,
strip.background.y = element_blank(),
axis.title.y = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "100 Mbp")
#
mp2 <- ggplot(x1, aes(x = `-log10(p)_Exp`, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_point(size = 0.5, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25, color = alpha("white", 0)) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_abline() +
scale_x_continuous(breaks = mybreaks, labels = mylabels) +
scale_y_continuous(breaks = mybreaks, labels = mylabels) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
facet_grid(Facet ~ "QQ", scales = "free_y") +
theme_AGL +
theme(legend.position = "none",
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
strip.text.y = ggtext::element_markdown()) +
labs(y = NULL, x = "Expected")
#
mp <- ggarrange(mp1, mp2, ncol = 2, widths = c(7,1.5), align = "h",
legend = "bottom", common.legend = T)
#
ggsave(paste0("Additional/CV/CV_", i, ".png"), mp, width = 12, height = 6)
}im1 <- magick::image_read("Additional/CV/CV_Ro17_DTF.png") %>%
magick::image_annotate("a)", size = 70, weight = 600, location = "+30+10")
im2 <- magick::image_read("Additional/CV/CV_Sp17_DTF.png") %>%
magick::image_annotate("b)", size = 70, weight = 600, location = "+30+10")
im <- magick::image_append(c(im1, im2), stack = T)
magick::image_write(im, "Figure_04.png")Figure 5
expttraits <- c("Su17_DTF", "Ba17_DTF", "Ne17_DTF", "It17_DTF")
expts <- substr(expttraits, 1, 4)
#
markers <- myMs[c(1,2)]
gg <- myG %>% filter(rs %in% markers) %>%
column_to_rownames("rs") %>%
select(11:ncol(.)) %>%
t() %>% as.data.frame()
for(i in 1:nrow(gg)) {
gg$Markers[i] <- paste(gg[i,1:length(markers)], collapse = " - ")
}
gg <- gg %>% rownames_to_column("Name")
#
yy <- myY %>% select(Name, expttraits) %>%
gather(ExptTrait, Value, 2:ncol(.)) %>%
separate(ExptTrait, c("Expt","Trait"), remove = F)
#
xx <- left_join(yy, gg, by = "Name") %>%
left_join(myLDP, by = "Name") %>%
mutate(Expt = factor(Expt, levels = expts))
#
x1 <- xx %>%
filter(!grepl("M|R|W|S|Y|K|N", Markers), ExptTrait %in% expttraits) %>%
arrange(Markers) %>%
mutate(Markers = factor(Markers, levels = c("C - A","C - C","T - A","T - C")))
#
mp1 <- ggplot(x1, aes(x = Markers, y = Value)) +
geom_quasirandom(aes(color = Cluster, key1 = Name, key2 = Origin)) +
facet_wrap(Expt ~ ., ncol = 6, scales = "free_y") +
scale_color_manual(values = myColors_Cluster) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3))) +
theme_AGL +
theme(legend.position = "bottom",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5)) +
labs(x = "Lcu.2RBY.Chr2p42556949 - ( C | T )\nLcu.2RBY.Chr5p1069654 - ( A | C )",
y = "DTF")
#
markers <- myMs[3]
gg <- myG %>% filter(rs %in% markers) %>%
column_to_rownames("rs") %>%
select(11:ncol(.)) %>%
t() %>% as.data.frame()
for(i in 1:nrow(gg)) {
gg$Markers[i] <- paste(gg[i,1:length(markers)], collapse = " - ")
}
gg <- gg %>% rownames_to_column("Name")
#
yy <- myY %>% select(Name, expttraits) %>%
gather(ExptTrait, Value, 2:ncol(.)) %>%
separate(ExptTrait, c("Expt","Trait"), remove = F)
#
xx <- left_join(yy, gg, by = "Name") %>%
left_join(myLDP, by = "Name") %>%
mutate(Expt = factor(Expt, levels = expts))
#
x1 <- xx %>%
filter(!grepl("M|R|W|S|Y|K|N", Markers)) %>%
arrange(Markers)
#
mp2 <- ggplot(x1 %>% filter(ExptTrait %in% expttraits), aes(x = Markers, y = Value)) +
geom_quasirandom(aes(color = Cluster, key1 = Name, key3 = Origin)) +
facet_wrap(Expt ~ ., ncol = 6, scales = "free_y") +
scale_color_manual(values = myColors_Cluster) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3))) +
theme_AGL +
theme(legend.position = "bottom",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5)) +
labs(x = "Lcu.2RBY.Chr6p2528817",
y = "DTF")
#
markers <- myMs[4]
gg <- myG %>% filter(rs %in% markers) %>%
column_to_rownames("rs") %>%
select(11:ncol(.)) %>%
t() %>% as.data.frame()
for(i in 1:nrow(gg)) {
gg$Markers[i] <- paste(gg[i,1:length(markers)], collapse = " - ")
}
gg <- gg %>% rownames_to_column("Name")
#
yy <- myY %>% select(Name, expttraits) %>%
gather(ExptTrait, Value, 2:ncol(.)) %>%
separate(ExptTrait, c("Expt","Trait"), remove = F)
#
xx <- left_join(yy, gg, by = "Name") %>%
left_join(myLDP, by = "Name") %>%
mutate(Expt = factor(Expt, levels = expts))
#
x1 <- xx %>%
filter(!grepl("M|R|W|S|Y|K|N", Markers)) %>%
arrange(Markers)
#
mp3 <- ggplot(x1 %>% filter(ExptTrait %in% expttraits), aes(x = Markers, y = Value)) +
geom_quasirandom(aes(color = Cluster, key1 = Name, key3 = Origin)) +
facet_wrap(Expt ~ ., ncol = 6, scales = "free_y") +
scale_color_manual(values = myColors_Cluster) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3))) +
theme_AGL +
theme(legend.position = "bottom",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5)) +
labs(x = "Lcu.2RBY.Chr6p307256203",
y = "DTF")
#
mp <- ggarrange(mp1, mp2, mp3, ncol = 1, labels = c("a)","b)","c)"),
heights = c(1.2,1,1), legend = "bottom", common.legend = T)
#
ggsave("Figure_05.png", mp, width = 8, height = 8)
ggsave("Additional/Figure_05_a.png", mp1, width = 8, height = 4)
ggsave("Additional/Figure_05_b.png", mp2, width = 8, height = 4)
ggsave("Additional/Figure_05_c.png", mp3, width = 8, height = 4)
#
mp1 <- plotly::ggplotly(mp1)
htmlwidgets::saveWidget(plotly::as_widget(mp1),
paste0("Additional/Figure_05_a_plotly.html"),
selfcontained = T)
#
mp2 <- plotly::ggplotly(mp2)
htmlwidgets::saveWidget(plotly::as_widget(mp2),
paste0("Additional/Figure_05_b_plotly.html"),
selfcontained = T)
#
mp3 <- plotly::ggplotly(mp3)
htmlwidgets::saveWidget(plotly::as_widget(mp3),
paste0("Additional/Figure_05_c_plotly.html"),
selfcontained = T)Supplemental Figure 1
gg_gwas_cv_summary <- function(filename, myTs, hlines, width, height) {
#
myTs <- c(myTs, paste0(myTs,"-b"), paste0(myTs, "-c"))
myLs <- gsub("-b","-*b*", myTs)
myLs <- gsub("-c","-*c*", myLs)
#
me <- c("<b style='color:black'>{ExptTrait}</b>",
"<b style='color:darkgreen'>{ExptTrait}</b>",
"<b style='color:darkorange3'>{ExptTrait}</b>",
"<b style='color:darkblue'>{ExptTrait}</b>")
#
e1 <- c("Ro16","Ro17","Su16","Su17","Su18","Us18")
e2 <- c("Mo16","Mo17","Sp16","Sp17","It16","It17")
e3 <- c("In16","In17","Ba16","Ba17","Ne16","Ne17")
myMacroEnvs <- c("Temperate", "South Asia", "Mediterranean")
#
xx <- myP %>%
filter(paste(Expt, Trait, sep = "_") %in% myTs) %>%
mutate(ExptTrait = paste(Expt, Trait, sep = "_"),
CV = "No Covariate",
CV = ifelse(grepl("-b", Trait), "CV = *b*", CV),
CV = ifelse(grepl("-c", Trait), "CV = *c*", CV),
CV = factor(CV, levels = c("CV = *c*","No Covariate","CV = *b*")),
Trait = gsub("-b|-c", "", Trait),
ExptTrait = gsub("-b|-c", "", ExptTrait),
Model = factor(Model, levels = myModels),
MacroEnv = NA,
MacroEnv = ifelse(Expt %in% e1, "Temperate", MacroEnv),
MacroEnv = ifelse(Expt %in% e2, "Mediterranean", MacroEnv),
MacroEnv = ifelse(Expt %in% e3, "South Asia", MacroEnv),
MacroEnv = ifelse(is.na(MacroEnv), "Multi Environment", MacroEnv),
cat_col = NA,
cat_col = ifelse(MacroEnv == "Multi Environment", glue::glue(me[1]), cat_col),
cat_col = ifelse(MacroEnv == "Temperate", glue::glue(me[2]), cat_col),
cat_col = ifelse(MacroEnv == "South Asia", glue::glue(me[3]), cat_col),
cat_col = ifelse(MacroEnv == "Mediterranean", glue::glue(me[4]), cat_col),
cat_col = factor(cat_col, levels = rev(cat_col[match(myTs, ExptTrait)])),
MacroEnv = factor(MacroEnv, levels = myMacroEnvs)) %>%
arrange(ExptTrait)
#
x1 <- xx %>% filter(`-log10(p)` > threshold)
x2 <- xx %>% filter(`-log10(p)` < threshold)
#
mp <- ggplot(x1, aes(x = Position / 100000000, y = cat_col,
shape = Model, fill = MacroEnv)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_point(data = x2, size = 1, color = "black", alpha = 0.5,
aes(key1 = SNP, key2 = FT_Genes, key3 = FT_Distances)) +
geom_point(size = 2, alpha = 0.5,
aes(key1 = SNP, key2 = FT_Genes, key3 = FT_Distances)) +
geom_hline(yintercept = hlines, alpha = 0.5) +
facet_grid(CV ~ Chromosome, drop = F, scales = "free_x", space = "free_x") +
scale_fill_manual(values = myColors_Macro[2:4], guide = "none") +
scale_shape_manual(values = c(22:25)) +
scale_y_discrete(drop = F) +
theme_AGL +
theme(legend.position = "bottom",
axis.text.y = element_markdown(),
strip.text.y = element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 4, fill = "black"))) +
labs(y = NULL, x = "100 Mbp")
#
ggsave(paste0(filename, ".png"), mp, width = width, height = height)
mp <- plotly::ggplotly(mp)
htmlwidgets::saveWidget(plotly::as_widget(mp),
paste0(filename, "_plotly.html"),
knitrOptions = list(fig.width = width, fig.height = height),
selfcontained = T)
}
gg_gwas_cv_summary(filename = "Supplemental_Figure_01",
myTs = c("Ro16_DTF", "Ro17_DTF", "Su16_DTF", "Su17_DTF", "Su18_DTF", "Us18_DTF",
"In16_DTF", "In17_DTF", "Ba16_DTF", "Ba17_DTF", "Ne16_DTF", "Ne17_DTF",
"Mo16_DTF", "Mo17_DTF", "Sp16_DTF", "Sp17_DTF", "It16_DTF", "It17_DTF"),
hlines = c(6.5,12.5), width = 12, height = 9)
gg_gwas_cv_summary(filename = "Additional/Additional_Figure_03",
myTs = c("Su17_Tf", "Ba17_Tf", "It17_Tf",
"Su17_Tb", "Ba17_Tb", "It17_Tb",
"Su17_Pf", "Ba17_Pf", "It17_Pf",
"Su17_Pc", "Ba17_Pc", "It17_Pc"),
hlines = c(3.5,6.5,9.6), width = 12, height = 9)
gg_gwas_cv_summary(filename = "Additional/Additional_Figure_04",
myTs = c("PCA_PC1", "PCA_PC2", "PCA_PC3",
"PTModel_a", "PTModel_b", "PTModel_c"),
hlines = c(3.5), width = 12, height = 6)Supplemental Figure 2
expttraits <- c("Su17_DTF", "Ba17_DTF", "Ne17_DTF",
"Sp17_DTF", "PCA_PC1", "PTModel_b")
chr <- 2
pos1 <- 35000000
pos2 <- 50000000
#
xx <- NULL
for(i in expttraits) {
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results/"))
fnames <- list.files("Results/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results/", j)) %>%
filter(Chromosome == chr, Position > pos1, Position < pos2) %>%
mutate(Model = mod, ExptTrait = i,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xx <- bind_rows(xx, xj)
}
}
i <- "Su17_DTF"
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results_c/"))
fnames <- list.files("Results_c/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results_c/", j)) %>%
filter(Chromosome == chr, Position > pos1, Position < pos2) %>%
mutate(Model = mod, ExptTrait = paste0(i,"-c"),
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xx <- bind_rows(xx, xj)
}
#
expttraits <- c("Su17_DTF", "Su17_DTF (CV=*c*)",
"Ba17_DTF", "Ne17_DTF", "Sp17_DTF", "PC1", "*b*")
xx <- xx %>%
mutate(ExptTrait = plyr::mapvalues(ExptTrait,
c("Su17_DTF-c", "PCA_PC1", "PTModel_b"),
c("Su17_DTF (CV=*c*)", "PC1", "*b*")),
ExptTrait = factor(ExptTrait, levels = expttraits),
Model = factor(Model, levels = myModels)) %>%
mutate(Chromosome = paste("Chromosome", Chromosome)) %>%
arrange(desc(Model))
#
x1 <- xx %>% filter(`-log10(p)` < threshold2)
x2 <- xx %>% filter(`-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(`-log10(p)` > threshold)
#
vi <- myGM %>% filter(Chromosome == 2) %>%
mutate(Chromosome = paste("Chromosome",Chromosome))
#
mp <- ggplot() +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_vline(data = vi, alpha = 0.5, color = "red",
aes(xintercept = Position / 1000000)) +
geom_point(data = x1, size = 0.75, color = alpha("white", 0),
aes(x = Position / 1000000, y = `-log10(p)`, shape = Model, fill = Model)) +
geom_point(data = x2, size = 1, color = alpha("white", 0),
aes(x = Position / 1000000, y = `-log10(p)`, shape = Model, fill = Model)) +
geom_point(data = x3, size = 1.5,
aes(x = Position / 1000000, y = `-log10(p)`, shape = Model, fill = Model)) +
facet_grid(ExptTrait ~ Chromosome, scales = "free", space = "free_x") +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(legend.position = "bottom",
legend.box = "vertical",
panel.grid = element_blank(),
axis.title.y = ggtext::element_markdown(),
strip.text.y = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "Mbp")
#
ggsave("Supplemental_Figure_02.png", mp, width = 5, height = 10)Supplemental Figure 3
expttraits <- c("PCA_PC3", "Su17_Tb", "It17_Tb", "Su17_Tf", "It17_Tf", "Mo16_DTF")
markers <- c("Lcu.2RBY.Chr6p307256203",
"Lcu.2RBY.Chr6p309410465",
"Lcu.2RBY.Chr6p306914970")
chr <- 6
pos1 <- 300000000
pos2 <- 315000000
colors <- c("darkgreen","darkred","darkblue","darkslategray4","khaki4")
#
xx <- NULL
for(i in expttraits) {
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results/"))
fnames <- list.files("Results/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results/", j)) %>%
filter(Chromosome == chr, Position > pos1, Position < pos2) %>%
mutate(Model = mod, ExptTrait = i,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xx <- bind_rows(xx, xj)
}
}
#
expttraits <- c("PC3", "Su17_*Tb*", "It17_*Tb*", "Su17_*Tf*", "It17_*Tf*", "Mo16_DTF")
xx <- xx %>%
mutate(Model = factor(Model, levels = myModels),
ExptTrait = plyr::mapvalues(ExptTrait,
c("PCA_PC3", "Su17_Tb", "It17_Tb", "Su17_Tf", "It17_Tf"),
c("PC3", "Su17_*Tb*", "It17_*Tb*", "Su17_*Tf*", "It17_*Tf*")),
ExptTrait = factor(ExptTrait, levels = expttraits)) %>%
mutate(Chromosome = paste("Chromosome",Chromosome)) %>%
arrange(desc(Model))
#
x1 <- xx %>% filter(`-log10(p)` < threshold2)
x2 <- xx %>% filter(`-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(`-log10(p)` > threshold)
#
myGenes <- c("LcLWD1", "LcFTa1", "LcFTc")
vg <- myFTGenes %>%
filter(Name %in% myGenes) %>%
mutate(Name = plyr::mapvalues(Name, myGenes, paste0("*",myGenes,"*")),
Name = factor(Name, levels = paste0("*",myGenes,"*")),
Chromosome = paste("Chromosome",Chromosome))
#
mp <- ggplot() +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_vline(data = vg, alpha = 0.5, color = "red",
aes(xintercept = Position / 1000000, lty = Name)) +
geom_point(data = x1, size = 0.75, color = alpha("white", 0),
aes(x = Position / 1000000, y = -log10(P.value),
shape = Model, fill = Model)) +
geom_point(data = x2, size = 1, color = alpha("white", 0),
aes(x = Position / 1000000, y = -log10(P.value),
shape = Model, fill = Model)) +
geom_point(data = x3, size = 1.5,
aes(x = Position / 1000000, y = -log10(P.value),
shape = Model, fill = Model)) +
facet_grid(ExptTrait ~ Chromosome, scales = "free", space = "free_x") +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
scale_linetype_manual(name = "Genes", values = c(4,1,3)) +
theme_AGL +
theme(legend.position = "bottom",
legend.box = "vertical",
panel.grid = element_blank(),
axis.title.y = ggtext::element_markdown(),
strip.text.y = ggtext::element_markdown(),
legend.text = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "Mbp")
ggsave("Supplemental_Figure_03.png", mp, width = 5, height = 10)Supplemental Figure 4
library(genetics)
dna <- data.frame(stringsAsFactors = F,
Symbol = c("A", "C", "G", "T", "U",
"R", "Y", "S", "W", "K", "M", "N"),
Value = c("A/A","C/C","G/G","T/T","U/U",
"A/G","C/T","G/C","A/T","G/T","A/C","N/N") )
xx <- myG[,c(-2,-5,-6,-7,-8,-9,-10,-11)]
for(i in 4:ncol(xx)) {
xx[xx[,i]=="N", i] <- NA
xx[,i] <- plyr::mapvalues(xx[,i], dna$Symbol, dna$Value)
}
# Num = 50
# x <- xx
LD_Decay <- function(x, folder = "Additional/LD/", Chr = 1, Num = 200) {
xc <- x %>% filter(chrom == Chr)
xr <- xc %>% column_to_rownames("rs")
xr <- xr[,c(-1,-2)]
rr <- round(runif(Num, 1, nrow(xr)))
while(sum(duplicated(rr))>0) {
ra <- round(runif(sum(duplicated(rr)), 1, nrow(xr)))
rr <- rr[!duplicated(rr)]
rr <- c(rr, ra)
}
rr <- rr[order(rr)]
xr <- xr[rr,]
#
xi <- xr %>% t() %>% as.data.frame()
myLD <- genetics::LD(genetics::makeGenotypes(xi))
save(myLD, file = paste0(folder, "LD_Chrom_", Chr, "_", Num, ".Rdata"))
}
#
LD_Decay(x = xx, Chr = 1, Num = 1000)
LD_Decay(x = xx, Chr = 2, Num = 1000)
LD_Decay(x = xx, Chr = 3, Num = 1000)
LD_Decay(x = xx, Chr = 4, Num = 1000)
LD_Decay(x = xx, Chr = 5, Num = 1000)
LD_Decay(x = xx, Chr = 6, Num = 1000)
LD_Decay(x = xx, Chr = 7, Num = 1000)# Create user-defined function
movingAverage <- function(x, n = 5) {
stats::filter(x, rep(1 / n, n), sides = 2)
}
#
xx <- NULL
for(i in 1:7) {
xc <- myG %>% filter(chrom == i)
load(paste0("Additional/LD/LD_Chrom_",i,"_1000.Rdata"))
xi <- myLD$`R^2` %>% as.data.frame() %>%
rownames_to_column("SNP1") %>%
gather(SNP2, LD, 2:ncol(.)) %>%
filter(!is.na(LD)) %>%
mutate(Chr = i,
SNP1_d = plyr::mapvalues(SNP1, xc$rs, xc$pos, warn_missing = F),
SNP2_d = plyr::mapvalues(SNP2, xc$rs, xc$pos, warn_missing = F),
Distance = as.numeric(SNP2_d) - as.numeric(SNP1_d)) %>%
arrange(Distance, rev(LD))
#
xii <- xi %>% filter(Distance < 1000000)
myloess <- stats::loess(LD ~ Distance, data = xii, span=0.50)
#
xi <- xi %>%
mutate(Moving_Avg = movingAverage(LD, n = 100),
Loess = ifelse(Distance < 1000000, predict(myloess), NA))
#
xx <- bind_rows(xx, xi)
}
#
x1 <- xx %>% group_by(Chr) %>%
summarise(Mean_LD = mean(Moving_Avg, na.rm = T))
x2 <- xx %>% filter(Loess < 0.2) %>% group_by(Chr) %>%
summarise(Threshold_0.2 = min(Distance, na.rm = T))
x3 <- xx %>% left_join(x1, by = "Chr") %>%
group_by(Chr) %>% filter(Moving_Avg < Mean_LD) %>%
summarise(Threshold_0.1 = min(Distance, na.rm = T))
myChr <- x1 %>%
left_join(x2, by = "Chr") %>%
left_join(x3, by = "Chr")
#
xx <- left_join(xx, myChr, by = "Chr") %>%
mutate(Chr = as.factor(Chr))
#
#mytitle <- paste(myChr$Chr, myChr$Threshold_0.2, sep = "-", collapse = "; ")
#
mp1 <- ggplot(xx, aes(x = Distance/1000000, y = Moving_Avg)) +
geom_line(size = 0.5, alpha = 0.5) +
geom_hline(data = myChr, aes(yintercept = Mean_LD), color = "red", lty = 2) +
geom_hline(yintercept = 0.2, color = "blue", lty = 2) +
scale_y_continuous(breaks = seq(0, 0.5, by = 0.1), limits = c(0,0.5)) +
facet_wrap(paste("Chr =", Chr) ~ ., ncol = 7, scales = "free_x") +
theme_AGL +
theme(legend.position = "none",
axis.title.y = ggtext::element_markdown()) +
labs(y = "R^2", x = "Mbp")
#
yy <- xx %>% filter(Distance < 1000000)
mp2 <- ggplot(yy, aes(x = Distance/1000, y = Loess)) +
geom_line(aes(y = Moving_Avg), alpha = 0.5) +
geom_line() +
geom_hline(data = myChr, aes(yintercept = Mean_LD), color = "red", lty = 2) +
geom_hline(yintercept = 0.2, color = "blue", lty = 2) +
scale_y_continuous(breaks = seq(0, 0.5, by = 0.1), limits = c(0,0.5)) +
facet_wrap(paste("Chr =", Chr) + paste("Threshold =", Threshold_0.2) ~ ., ncol = 7, scales = "free_x") +
geom_vline(data = myChr, lty = 2, size = 0.3, alpha = 0.8,
aes(xintercept = Threshold_0.2/1000)) +
theme_AGL +
theme(legend.position = "none",
axis.title.y = ggtext::element_markdown()) +
labs(y = "R^2", x = "Kbp")
#
mp <- ggarrange(mp1, mp2, ncol = 1, labels = c("a)", "b)"))
#
ggsave(paste0("Supplemental_Figure_04.png"), mp , width = 12, height = 8)Supplemental Figure 5
Supplemental Table 2
Additional Plots
Phenotype Data
xx <- myY %>%
gather(ExptTrait, Value, 2:ncol(.)) %>%
mutate(ExptTrait = factor(ExptTrait, levels = colnames(myY)[-1]))
#
mp <- ggplot(xx, aes(x = Value)) +
geom_histogram(fill = "darkgreen", color = "black", alpha = 0.7) +
facet_wrap(ExptTrait~., ncol = 9, scales = "free") +
theme_AGL
ggsave("Additional/Additional_Figure_00.png", mp, width = 14, height = 7)Grouped Manhattan Plots
gg_Custom_Man <- function(expttraits = c("Su17_DTF", "Ba17_DTF", "Ne17_DTF", "It17_DTF"),
folder = "Results/", filename = "Fig01.png", title = filename,
markers = myMs, height = 10) {
#
#vt <- myP %>% filter(paste(Expt, Trait, sep = "_") %in% expttraits) %>%
# arrange(P.value) %>%
# filter(!duplicated(SNP)) %>%
# mutate(FT_Genes = ifelse(FT_Distance < 1000000, FT_Genes, NA),
# ExptTrait = paste(Expt, Trait, sep="_"))
xx <- NULL
xs <- NULL
for(i in expttraits) {
fnames <- grep(paste0(i,".GWAS.Results"), list.files(folder))
fnames <- list.files(folder)[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xsi <- GWAS_PeakTable(folder = folder, file = j) %>%
arrange(P.value) %>%
filter(!duplicated(paste(FT_Genes, Model))) %>%
mutate(FT_Genes = gsub(" ; ", "\n", FT_Genes))
xj <- read.csv(paste0(folder, j)) %>%
mutate(Model = mod, ExptTrait = i,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xx <- bind_rows(xx, xj)
if(nrow(xsi)>0) { xs <- bind_rows(xs, xsi) }
}
}
#
xx <- xx %>%
mutate(Model = factor(Model, levels = myModels),
ExptTrait = factor(ExptTrait, levels = expttraits)) %>%
arrange(desc(Model))
#
xx <- xx %>% mutate(`-log10(p)` = ifelse(`-log10(p)` > 10, 10, `-log10(p)`))
xs <- xs %>% mutate(`-log10(p)` = ifelse(`-log10(p)` > 10, 10, `-log10(p)`),
ExptTrait = paste(Expt, Trait, sep = "_"))
#
x1 <- xx %>% filter(`-log10(p)` < threshold2)
x2 <- xx %>% filter(`-log10(p)` > threshold2, `-log10(p)` < threshold)
x3 <- xx %>% filter(`-log10(p)` > threshold)
#
mp <- ggplot(x1, aes(x = Position / 100000000, y = `-log10(p)`,
shape = Model, fill = Model)) +
geom_vline(data = myGM, color = "red", alpha = 0.6,
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
ggrepel::geom_text_repel(data = xs, aes(label = FT_Genes, ggroup = ExptTrait),
size = 1.5, nudge_y = 0.5) +
facet_grid(ExptTrait ~ Chromosome, scales = "free_x", space = "free_x") +
scale_x_continuous(breaks = 1:7) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(legend.position = "bottom") +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(title = title, y = NULL, x = "100 Mbp")
#
ggsave(paste0("Additional/Man_Grouped/", filename, ".png"),
mp, width = 10, height = height)
}
# DTF
gg_Custom_Man(expttraits = c("Ro16_DTF", "Ro17_DTF", "Su16_DTF",
"Su17_DTF", "Su18_DTF", "Us18_DTF"),
filename = "Man_DTF_Temperate")
gg_Custom_Man(expttraits = c("It16_DTF", "It17_DTF", "Sp16_DTF",
"Sp17_DTF", "Mo16_DTF", "Mo17_DTF"),
filename = "Man_DTF_Mediterranean")
gg_Custom_Man(expttraits = c("Ne16_DTF", "Ne17_DTF", "Ba16_DTF",
"Ba17_DTF", "In16_DTF", "In17_DTF"),
filename = "Man_DTF_SouthAsia")
# Tf & Tb
gg_Custom_Man(expttraits = c("Su17_Tf", "Ba17_Tf", "It17_Tf",
"Su17_Tb", "Ba17_Tb", "It17_Tb"),
filename = "Man_Tf_Tb")
# Pf & Pc
gg_Custom_Man(expttraits = c("Su17_Pf", "Ba17_Pf", "It17_Pf",
"Su17_Pc", "Ba17_Pc", "It17_Pc"),
filename = "Man_Pf_Pc")
# PCA & abc
gg_Custom_Man(expttraits = c("PCA_PC1", "PCA_PC2", "PCA_PC3",
"PTModel_a", "PTModel_b", "PTModel_c"),
filename = "Man_PCA_abc")
# DTF - c
gg_Custom_Man(expttraits = c("Ro16_DTF", "Ro17_DTF", "Su16_DTF",
"Su17_DTF", "Su18_DTF", "Us18_DTF"),
folder = "Results_c/",
filename = "Man_DTF_c_Temperate")
gg_Custom_Man(expttraits = c("It16_DTF", "It17_DTF", "Sp16_DTF",
"Sp17_DTF", "Mo16_DTF", "Mo17_DTF"),
folder = "Results_c/",
filename = "Man_DTF_c_Mediterranean")
gg_Custom_Man(expttraits = c("Ne16_DTF", "Ne17_DTF", "Ba16_DTF",
"Ba17_DTF", "In16_DTF", "In17_DTF"),
folder = "Results_c/",
filename = "Man_DTF_c_SouthAsia")
# DTF - b
gg_Custom_Man(expttraits = c("Ro16_DTF", "Ro17_DTF", "Su16_DTF",
"Su17_DTF", "Su18_DTF", "Us18_DTF"),
folder = "Results_b/",
filename = "Man_DTF_b_Temperate")
gg_Custom_Man(expttraits = c("It16_DTF", "It17_DTF", "Sp16_DTF",
"Sp17_DTF", "Mo16_DTF", "Mo17_DTF"),
folder = "Results_b/",
filename = "Man_DTF_b_Mediterranean")
gg_Custom_Man(expttraits = c("Ne16_DTF", "Ne17_DTF", "Ba16_DTF",
"Ba17_DTF", "In16_DTF", "In17_DTF"),
folder = "Results_b/",
filename = "Man_DTF_b_SouthAsia")Facetted & Multi-Modeled Manhattan Plots
gg_Man <- function(folder = "Results/", expttrait = "Ro17_DTF", suffix = NULL,
colors = c("darkgreen","darkgoldenrod3","darkgreen","darkgoldenrod3",
"darkgreen", "darkgoldenrod3","darkgreen")) {
#
expt <- substr(expttrait, 1, gregexpr("_", expttrait)[[1]][1]-1 )
trait <- substr(expttrait, gregexpr("_", expttrait)[[1]][1]+1, nchar(expttrait) )
fnames <- grep(paste0(expttrait, ".GWAS.Results"), list.files(folder))
fnames <- list.files(folder)[fnames]
xx <- NULL
xs <- NULL
for(i in fnames) {
trt <- substr(i, gregexpr("GAPIT.", i)[[1]][1]+6,
gregexpr(".GWAS.Results.csv", i)[[1]][1]-1 )
mod <- substr(trt, 1, gregexpr("\\.", trt)[[1]][1]-1 )
trt <- substr(trt, gregexpr("\\.", trt)[[1]][1]+1, nchar(trt) )
#
xsi <- GWAS_PeakTable(folder = folder, file = i) %>%
arrange(P.value) %>%
filter(!duplicated(paste(FT_Genes, Model))) %>%
mutate(FT_Genes = gsub(" ; ", "\n", FT_Genes))
xi <- read.csv(paste0(folder, i))
#
if(sum(colnames(xi)=="nobs")>0) { xi <- select(xi, -nobs) }
xi <- xi %>%
mutate(Model = mod,
`-log10(p)` = -log10(P.value),
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)),
`-log10(p)_FDR` = -log10(FDR_Adjusted_P.values))
xx <- bind_rows(xx, xi)
if(nrow(xsi)>0) { xs <- bind_rows(xs, xsi) }
}
#
xs2 <- xs %>% select(SNP, Model, FT_Genes)
xx <- xx %>%
mutate(Sig = ifelse(-log10(P.value) > threshold2, "Suggested", "Not Significant"),
Sig = ifelse(-log10(P.value) > threshold, "Significant", Sig)) %>%
left_join(xs2, by = c("SNP", "Model")) %>%
mutate(Model = factor(Model, levels = myModels)) %>%
arrange(desc(Model))
#
xx <- xx %>% mutate(`-log10(p)` = ifelse(`-log10(p)` > 10, 10, `-log10(p)`))
#
x1 <- xx %>% filter(-log10(P.value) > 0.5, -log10(P.value) < threshold2)
x2 <- xx %>% filter(-log10(P.value) > threshold2 & -log10(P.value) < threshold)
x3 <- xx %>% filter(-log10(P.value) > threshold)
# Man plot
mp1 <- ggplot(xx, aes(x = Position / 100000000, y = -log10(P.value))) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_point(aes(color = factor(Chromosome)), pch = 16, size = 0.5) +
geom_point(data = x2, pch = 18, size = 0.75, color = "darkred") +
geom_point(data = x3, pch = 23, size = 1.25, color = "black", fill = "darkred") +
ggrepel::geom_text_repel(aes(label = FT_Genes), size = 1.5, nudge_y = 0.5) +
facet_grid(Model ~ Chromosome, scales = "free") +
scale_x_continuous(breaks = 0:7) +
scale_color_manual(values = colors) +
theme_AGL +
theme(legend.position = "none",
axis.title.y = ggtext::element_markdown()) +
labs(title = paste0(expttrait, suffix),
y = "-log<sub>10</sub>(*p*)", x = "100 Mbp")
# QQ plot
mp2 <- ggplot(x1, aes(y = `-log10(p)`, x = `-log10(p)_Exp`)) +
geom_point(pch = 16, size = 0.75, color = colors[1]) +
geom_point(data = x2, pch = 18, size = 1, color = "darkred") +
geom_point(data = x3, pch = 23, size = 1.25, color = "black", fill = "darkred") +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_abline() +
facet_grid(Model ~ "QQ", scales = "free_y") +
theme_AGL +
labs(title = "", y = NULL, x = "Expected")
# Append and save plots
mp <- ggarrange(mp1, mp2, ncol = 2, widths = c(4,1), align = "h")
ggsave(paste0("Additional/Man_Facet/ManQQ_", expt, "_", trait, suffix, ".png"),
mp, width = 10, height = 7)
#
xs <- xs %>% arrange(P.value) %>% filter(!duplicated(FT_Genes))
# Phenotype Plots
mp1.1 <- ggplot(myY %>% select(Name, Value=expttrait)) +
geom_histogram(aes(x = Value), fill = colors[1], alpha = 0.8) +
facet_grid(. ~ paste0(expttrait, suffix)) +
theme_AGL +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
labs(x = "Value", y = "Count")
mp1.2 <- ggplot(myY %>% select(Name, Value=expttrait)) +
stat_ecdf(aes(x = Value), color = colors[1], alpha = 0.8, lwd = 1.5) +
facet_grid(. ~ "Accum. Dist.") +
theme_AGL +
labs(x = "Density", y = "Accumulation")
# Manhattan Plot
mp2 <- ggplot(x1, aes(x = Position / 100000000, y = -log10(P.value),
shape = Model, fill = Model)) +
geom_vline(data = myGM, alpha = 0.5, color = "red",
aes(xintercept = Position / 100000000)) +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
ggrepel::geom_text_repel(data = xs, aes(label = FT_Genes), size = 1.5, nudge_y = 0.5) +
facet_grid(. ~ Chromosome, scales = "free") +
scale_x_continuous(breaks = 0:7) +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(axis.title.y = ggtext::element_markdown()) +
guides(shape = guide_legend(override.aes = list(size = 3))) +
labs(y = "-log<sub>10</sub>(*p*)", x = "Mbp")
# QQ plot
mp3.1 <- ggplot(x1, aes(y = `-log10(p)`, x = `-log10(p)_Exp`, shape = Model, fill = Model)) +
geom_point(size = 0.5, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_abline() +
scale_fill_manual(name = NULL, values = myColors_Model) +
scale_shape_manual(name = NULL, values = 22:25) +
facet_grid(. ~ "QQ", scales = "free_y") +
theme_AGL +
theme(legend.position = "none") +
labs(title = "", y = NULL, x = "Expected")
# MAF plot
mp3.2 <- ggplot(x1, aes(y = `-log10(p)`, x = maf * 100, shape = Model, fill = Model)) +
geom_point(size = 0.3, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
geom_hline(yintercept = threshold, color = "red") +
geom_hline(yintercept = threshold2, color = "blue") +
geom_point(size = 0.3, color = alpha("white", 0)) +
facet_grid(. ~ "MAF", scales = "free_y") +
scale_fill_manual(name = NULL, values = myColors_Model) +
scale_shape_manual(name = NULL, values = 22:25) +
theme_AGL +
theme(legend.position = "none") +
labs(y = NULL, x = "%")
# Append and save plots
mp1 <- ggarrange(mp1.1, mp1.2, ncol = 1, align = "v")
mp3 <- ggarrange(mp3.1, mp3.2, ncol = 1, align = "v")
mp <- ggarrange(mp1, mp2, mp3, ncol = 3, align = "h",
widths = c(1.5,7,1.5), legend = "top")
#
ggsave(paste0("Additional/Man_Multi/ManQQ_", expt, "_", trait, suffix, ".png"),
mp, width = 12, height = 4)
}
#
expttraits <- colnames(myY)[2:ncol(myY)]
for(i in expttraits) {
gg_Man(folder = "Results/", expttrait = i)
}
for(i in expttraits[!grepl("PTModel_b", expttraits)]) {
gg_Man(folder = "Results_b/", expttrait = i, suffix = "-b")
}
for(i in expttraits[!grepl("PTModel_c", expttraits)]) {
gg_Man(folder = "Results_c/", expttrait = i, suffix = "-c")
}Marker Plots
gg_Marker <- function(marker = myMs[1],
expttraits = c("Su17_DTF", "Ba17_DTF", "Ne17_DTF", "It17_DTF",
"Su17_Tf", "Su17_Pf", "PTModel_b", "PTModel_c")) {
#
gg <- myG %>% filter(rs == marker) %>%
column_to_rownames("rs") %>%
select(11:ncol(.)) %>%
t() %>% as.data.frame() %>%
rownames_to_column("Name") %>%
rename(Marker=2)
#
yy <- myY %>% select(Name, expttraits) %>%
gather(ExptTrait, Value, 2:ncol(.))
#
xx <- left_join(yy, gg, by = "Name") %>%
left_join(myLDP, by = "Name") %>%
mutate(ExptTrait = factor(ExptTrait, levels = expttraits))
#
x1 <- xx %>%
filter(!grepl("M|R|W|S|Y|K|N", Marker), ExptTrait %in% expttraits) %>%
arrange(Marker)
#
fg <- myP %>% filter(SNP == marker) %>% slice(1)
fg <- paste(fg$FT_Genes, fg$FT_Distances, sep = "\n")
mp <- ggplot(x1, aes(x = Marker, y = Value)) +
geom_quasirandom(aes(color = Cluster)) +
facet_wrap(ExptTrait ~ ., ncol = 4, scales = "free_y") +
scale_color_manual(values = myColors_Cluster) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 3))) +
theme_AGL +
theme(legend.position = "bottom",
panel.grid.major.x = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5)) +
labs(title = marker, caption = fg,
x = NULL, y = "DTF")
mp
}
#
i<-myMs[1]
for(i in myMs) {
mp <- gg_Marker(marker = i)
ggsave(paste0("Additional/Markers/Markers_",i,".png"), mp, width = 8, height = 6)
}
pdf("Additional/Markers/Markers_Chr_1.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 1)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()
#
pdf("Additional/Markers/Markers_Chr_2.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 2)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()
#
pdf("Additional/Markers/Markers_Chr_3.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 3)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()
#
pdf("Additional/Markers/Markers_Chr_4.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 4)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()
#
pdf("Additional/Markers/Markers_Chr_5.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 5)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()
#
pdf("Additional/Markers/Markers_Chr_6.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 6)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()
#
pdf("Additional/Markers/Markers_Chr_7.pdf", width = 8, height = 6)
myPi <- myP %>% filter(Chromosome == 7)
for(i in unique(myPi$SNP)) { print(gg_Marker(marker = i)) }
dev.off()Flowering Time Genes GWAS Results
pdfGenePlots <- function(GeneName = "LcELF3a") {
myGene <- myFTGenes %>% filter(Name == GeneName)
expttraits <- c(
"PCA_PC1", "PCA_PC2", "PCA_PC3",
"PTModel_a", "PTModel_b", "PTModel_c",
"Ro16_DTF", "Ro17_DTF", "Su16_DTF", "Su17_DTF", "Su18_DTF", "Us18_DTF",
"In16_DTF", "In17_DTF", "Ba16_DTF", "Ba17_DTF", "Ne16_DTF", "Ne17_DTF",
"Mo16_DTF", "Mo17_DTF", "Sp16_DTF", "Sp17_DTF", "It16_DTF", "It17_DTF",
"Su17_Tf", "Ba17_Tf", "It17_Tf",
"Su17_Tb", "Ba17_Tb", "It17_Tb",
"Su17_Pf", "Ba17_Pf", "It17_Pf",
"Su17_Pc", "Ba17_Pc", "It17_Pc")
#
pdf(paste0("Additional/",myGene$Name, ".pdf"), width = 10, height = 6)
for(i in expttraits) {
x1 <- NULL; xb <- NULL; xc <- NULL
#
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results/"))
fnames <- list.files("Results/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results/", j)) %>%
mutate(Model = mod, Facet = i,
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
x1 <- bind_rows(x1, xj)
}
#
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results_b/"))
fnames <- list.files("Results_b/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results_b/", j)) %>%
mutate(Model = mod, Facet = "CV = *b*",
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xb <- bind_rows(xb, xj)
}
#
fnames <- grep(paste0(i,".GWAS.Results"), list.files("Results_c/"))
fnames <- list.files("Results_c/")[fnames]
for(j in fnames) {
mod <- substr(j, gregexpr("\\.", j)[[1]][1]+1, nchar(j))
mod <- substr(mod, 1, gregexpr("\\.", mod)[[1]][1]-1 )
xj <- read.csv(paste0("Results_c/", j)) %>%
mutate(Model = mod, Facet = "CV = *c*",
`-log10(p)_Exp` = -log10((rank(P.value, ties.method="first")-.5)/nrow(.)))
xc <- bind_rows(xc, xj)
}
#
xx <- bind_rows(x1, xb, xc) %>%
filter(Chromosome == myGene$Chromosome,
Position > myGene$Start - 10000000,
Position < myGene$Start + 10000000) %>%
mutate(Model = factor(Model, levels = myModels),
Facet = factor(Facet, levels = c("CV = *c*", i, "CV = *b*"))) %>%
arrange(Facet, rev(Model))
#
x1 <- xx %>% filter(-log10(P.value) < threshold2)
x2 <- xx %>% filter(-log10(P.value) > threshold2, -log10(P.value) < threshold)
x3 <- xx %>% filter(-log10(P.value) > threshold)
#
print(ggplot(x1, aes(x = Position / 100000000, y = -log10(P.value),
shape = Model, fill = Model)) +
geom_vline(data = myGene, xintercept = myGene$Start / 100000000) +
geom_hline(yintercept = threshold, color = "red", alpha = 0.7) +
geom_hline(yintercept = threshold2, color = "blue", alpha = 0.7) +
geom_point(size = 0.5, color = alpha("white", 0)) +
geom_point(data = x2, size = 0.75, color = alpha("white", 0)) +
geom_point(data = x3, size = 1.25) +
facet_grid(Facet ~ Chromosome, scales = "free", space = "free_x") +
scale_fill_manual(values = myColors_Model) +
scale_shape_manual(values = 22:25) +
theme_AGL +
theme(axis.title.y = ggtext::element_markdown(),
strip.text.y = ggtext::element_markdown()) +
labs(title = i, subtitle = GeneName,
y = "-log<sub>10</sub>(*p*)", x = "100 Mbp")
)
}
dev.off()
}
pdfGenePlots(GeneName = "LcELF3a")
pdfGenePlots(GeneName = "LcFTb1")
pdfGenePlots(GeneName = "LcFTa1")
pdfGenePlots(GeneName = "LcGI")© Derek Michael Wright